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Abstract 

This article corrsiders the detection of point sources in two dimensional astro¬ 
nomical images. The detection scheme we propose is based on peak statistics. We 
discuss the example of the detection of far galaxies in Cosmic Microwave Background 
experiments throughout the paper, although the method we present is totally gen¬ 
eral and can be used in many other fields of data analysis. We assume sources with 
a Gaussian profile -that is a fair approximation of the profile of a point source con¬ 
volved with the detector beam in microwave experiments- on a background modeled 
by a homogeneous and isotropic Gaussian random field characterized by a scale-free 
power spectrum. Point sources are enhanced with respect to the background by 
means of linear filters. After filtering, we identify local maxima and apply our de¬ 
tection scheme, a Neyman-Pearson detector that defines our region of acceptance 
based on the a priori pdf of the sources and the ratio of number densities. We study 
the different performances of some linear filters that have been used in this context 
in the literature: the Mexican Hat wavelet, the matched filter and the scale-adaptive 
filter. We consider as well an extension to two dimensions of the biparametric scale 
adaptive filter (BSAF). The BSAF depends on two parameters which are deter¬ 
mined by maximizing the number density of real detections while fixing the number 
density of spurious detections. For our detection criterion the BSAF outperforms 
the other filters in the interesting case of white noise. 

keywords: methods: analytical - methods: data analysis - techniques: image 
processing 


1 


1 INTRODUCTION 

A very challenging aspect of data analysis in astronomy is the detection of pointlike 
sources embedded in one and two dimensional images. Some common examples are 
the separation of individual stars in crowded optical images, the identihcation of 
emission and absorption lines in noisy one dimensional spectra and the detection of 
faint extragalactic objects at microwave frequencies. This latter case, for example, 
is one of the most critical issues for the new generation of experiments that observe 
the Cosmic Microwave Background (CMB). 

The CMB is the remnant of the radiation that filled the Universe immediately 
after the Big Bang. This weak radiation can provide us with answers to one of 
the most important set of questions asked in modern science - how the Universe 
did begin, how it evolved to the state we observe today, and how it will continue 
to evolve in the future. Unfortunately, we do not measure the CMB alone but a 
mixture of it with instrumental noise and other astrophysical radiations that are 
usually referred to as foregrounds. 

Some foregrounds are due to our own Galaxy, for example the thermal emission 
due to dust grains in the Galactic plane or the synchrotron emission by relativis¬ 
tic electrons moving along the Galactic magnetic field. These foregrounds appear 
as diffuse emission in the sky, and their spectral behaviours (the way the emission 
scales from one wavelength of observation to another) are reasonably well known. 
Another foreground with a well known spectral behaviour is the Sunyaev-Zel’dovich 
effect, which is due to the hot gas contained in galaxy clusters that distorts the en¬ 
ergy distribution of CMB photons. Foreground emissions carry information about 
the Galaxy structure, composition and physical parameters as well as about the 
number, distribution and evolution of galaxy clusters that map the distribution 
of matter in the Universe. Therefore, the study of the different foregrounds has 
great scientific relevance by itself. In order to properly study the CMB and the 
different foregrounds it is mandatory to separate the signals (components) that are 
mixed in the observations. This can be done by observing the sky at a number of 
frequencies at least as big as the number of components and then applying some sta¬ 
tistical component separation method in order to recover the different astrophysical 
signals. Several component separation techniques have been suggested, including 
blind (Baccigalupi et al. 2000, Maino et al. 2002, Delabrouille et al. 2003), semi¬ 
blind (Bedini et al. 2004) and non-blind (Hobson et al. 1998, Bouchet and Gispert 
1999, Stolyarov et al. 2002, Barreiro et al. 2004) approaches. 

Another important foreground is due to the emission of far galaxies. Since the 
typical angular size of the galaxies in the sky is a few arcseconds and the angular 
resolution of the microwave detectors is typically greater than a few arcminutes^, 
galaxies appear as points to the detector, which is unable to resolve their inner 
structure. Therefore, they are usually referred to as extragalactic point sources 

^For example, the upcoming ESA’s Planck satellite will have angular resolutions ranging from 5 
arcminutes (for the 217-857 GHz channels) to 33 arcminutes (for the 30 GHz channel). 
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(EPS) in the CMB jargon. Note that, however, they do not appear as points in the 
images but as the convolution of a pointlike impulse with the angular response of the 
detector (“beam”). The instruments (radiometers and bolometers) that are used 
in CMB experiments have angular responses that are approximately Gaussian and 
therefore EPS appear as small Gaussian (or nearly Gaussian) spots in the images^. 

The problem with EPS is that galaxies are a very heterogeneous bundle of ob¬ 
jects, from the radio galaxies that emit most of their radiation in the low frequency 
part of the electromagnetic spectrum to the dusty galaxies that emit mainly in the 
infrared (Toffolatti et al. 1998, Guiderdoni et al. 1998, Tucci et al. 2004). This 
makes it impossible to consider all of them as a single foreground to be separated 
from the other by means of multi-wavelength observations and statistical component 
separation techniques. EPS constitute an important contaminant in CMB studies 
at small angular scales (Toffolatti et al. 1998), affecting the determination of the 
CMB angular power spectrum and hampering the statistical study (e.g. the study 
of Gaussianity) of CMB and other foregrounds at such scales. Moreover, while there 
are good galaxy surveys at radio and infra-red frequencies, the microwave window 
of the electromagnetic spectrum is a practically unknown zone for extragalactic as¬ 
tronomy. Therefore, it is important to have detection techniques that are able to 
detect EPS with fluxes as low as possible. 

One possibility is to consider the EPS emission at each frequency as an addi¬ 
tional noise to be considered in the equations of a statistical component separation 
method. Once the algorithm has separated the different components, the resid¬ 
ual that is obtained by subtracting the output foregrounds from the original data 
should contain the EPS plus the instrumental noise and some amount of foreground 
residuals that remain due to a non-perfect separation. As an example, Eigure 1 
shows the residual at 30 GHz after applying a Maximum Entropy component sepa¬ 
ration algorithm (Hobson et al. 1999) to a 12.8 x 12.8 sqr deg simulated sky patch 
as would be observed by the Planck satellite. The brightest point sources can be 
clearly observed over the residual noise. However, fainter point sources are still 
masked by a residual noise that is approximately Gaussian and must be detected 
somehow. Besides, the situation is more complex because the presence of bright 
EPS in the data affects the performance of the component separation algorithms 
so the recovered components are contaminated by point sources in a way that is 
difficult to control. Therefore, any satisfactory method should detect and extract 
at least the bright sources before the component separation. Then, after separation 
some additional low intensity EPS could be detected from the residual maps such 
as the one in Eigure 1. 

Several techniques based on linear filters have been proposed in the literature for 
the detection of point sources in GMB data. Linear filtering techniques are suitable 

^It is also common to speak of compact sources, describing a source that is comparable to the size 
of the beam being used. Non-pointlike sources (such as large galaxy clusters with arcminute angular 
scales) will have more complicated responses when convolved with a beam, but if the source prohle is 
known it is always possible to apply the methods presented in this work. 
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Figure 1: Residual map of a 12.8 x 12.8 sqr deg sky patch at 30 GHz after the appli¬ 
cation of a Maximum Entropy component separation. The residual map is obtained 
by subtracting from the 30 GHz map the different components (CMB and foregrounds) 
given by the Maximum Entropy algorithm. Bright point sources appear as spots in the 
images whereas faint point sources are masked by the residual noise. 


4 


100 0 100 200 














for this problem because they can isolate structures with a given characteristic 
scale, as is the case of pointlike sources, while cancelling the contribution of diffuse 
foregrounds. Among the methods proposed in the literature, we emphasize the 
Mexican Hat Wavelet (MHW, Cayon et al. 2000, Vielva et al. 2001a,b,2003), the 
classic matched filter (MF, Tegmark and de Oliveira-Costa 1998), the Adaptive Top 
Hat Filter (Chiang et al. 2002) and the scale-adaptive filter (SAF, Sanz et al. 2001, 
Herranz et al. 2002). Moreover, linear filters can be used in combination with 
statistical component separation techniques in order to produce a more accurate 
separation of the different foregrounds (Vielva et al. 2001b). 

The goal of filtering is to enhance the contrast between the source to be detected 
and the background that masks it. For example, if we filter the image in Figure 1, 
assuming that the background can be characterised by a white noise, with the well 
known matched filter (see section 4.1) at the scale of the 30 GHz detector beam 
(FWHM=33 arcminutes) the signal to noise ratio of the sources increases by more 
than 25%. Therefore, a source whose signal to noise ratio was ~ 3 before filtering 
becomes a source with signal to noise ratio ~ 4 and will be easier to detect. 

After filtering, a deteetion rule is applied to the data in order to decide whether 
the source is present or not. The usual detection approach in astronomy is thresh¬ 
olding: for any given candidate (for example a local peak in the data) a positive 
detection is considered if the candidate has a signal to noise ratio greater than a 
certain threshold (in many astronomical applications, a typical value of this thresh¬ 
old is 5cr). This naive approach works fine for bright sources, but weak sources can 
be easily missed. 

More sophisticated detection schemes can use additional information in order to 
improve the detection. If the detection is performed by means of the study of the 
statistics of maxima in the images, such information includes not only the ampli¬ 
tude of the maxima but also spatial information related with the source profile, for 
example the derivatives of the intensity; in our approach we will consider the ampli¬ 
tude, the curvature and the shear of the sources (the last two quantities are given 
by the properties of the beam in the case of point sources) to discriminate between 
maxima of the background and real sources. Moreover, in some cases a priori in¬ 
formation on the distribution of intensity of the sources is known. We will therefore 
use a Neyman-Pearson detector that uses the three above mentioned elements of 
information (amplitude, curvature and shear) of the maxima as well as the a priori 
probability distribution of the sources. This technique has been successfully tested 
in images of one-dimensional fields (Lopez-Caniego et al. 2004a,b). In this work we 
will generalise it to two dimensions. 

The overview of this work is as follows: in section 2 we describe the statistics of 
the peaks for a two-dimensional Gaussian background in the absence and presence 
of a source. In section 3 we introduce the detection problem, define the region of 
acceptance and derive our detector. In section 4 we briefly review some of the linear 
filters proposed in the literature. In section 5 we describe a probability distribution 
of sources that is of interest and compare the performance of the filters, regarding 
our choice of detector. Finally, in section 6 we summarise our results. 
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2 PEAK STATISTICS 


In this section we will study the statistics of peaks for a two-dimensional Gaus¬ 
sian background in both the absence and presence of a source. We will focus on 
three quantities that define the properties of the peaks: the intensity of the field, 
the curvature and the shear at the position of the peak. The first quantity gives 
the amplitude of the peak. The curvature and the shear give information about 
the spatial structure of the peak and are related to its sharpness and eccentricity, 
respectively. 


2.1 Background 

Let us assume a two-dimensional (2D) background represented by a Gaussian ran¬ 
dom field ^(x) with average value {^{x)) = 0 and power spectrum P(g), 

{aQ)C{Q')) = P{q)Sh{Q-Q'), q=\Q\, (i) 


where ^{Q) is the Fourier transform of ^{x)^ and is the Dirac distribution in 2D. 

We are interested in the distribution of maxima of the background with respect to 
the three variables already mentioned: intensity, curvature and shear. Let us define 
the normalized field intensity v, the normalized curvature n and the normalized 
shear e as: 


n=-^. 

cro' 


K = 


Ai -l- A2 


e = 


0-2 


Ai — A2 

2(72 


( 2 ) 


where v G (— 00 , 00 ), k G [0, 00 ), e G [0,«:/2), Ai and A 2 are the eigenvalues of the 
negative Hessian matrix, and the (j„ are defined as 


1 

^ - I dq gi+2”P(<?). (3) 

The moment (Tq is equal to the dispersion of the field. 

The distribution of maxima of the background in one dimension (ID) with re¬ 
spect to the intensity and curvature (the shear is not defined in ID) was studied 
by Rice (1954). If we generalize it to 2D, including the shear, the expected number 
density of maxima per intervals (x, x + dx), {v, v -G dv), {k, k + dn) and (e, e -I- de) 
is given by 


nb{v, K, e) 


8\/3nb 


e{K^ 


4e^)e 


-4e^- 


2 


2(l-p^) ^ 


(4) 


^Throughout this paper we will use the following notation for the Fourier transform: the same symbol 
will be used for the real space and the Fourier space versions of a given function. The argument of the 
function will specify in each case which is the space we are referring to. For instance, f{q) will be the 
Fourier transform of the function f{x). 
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where hb is the the expected total number density of maxima 
maxima per unit area dx) 


rib = 


1 

47rv^ ’ 


(i.e. 


number of 


(5) 


and p and 9m are defined as 


Or, 


^2 croCr2 Oc <7i 


( 6 ) 


In the previous equations 9c and 9m are the coherence scale of the field and maxima, 
respectively. The formula in equation (4) can be derived from previous works (Bond 
and Efstathiou 1987, Barreiro et al. 1997). 


2.2 Background plus point source 

To the previous 2D background we add a source with a known spatial profile t{x) 
and an amplitude A, so that the intensity due to the source at a given position xq 
is ^s(al) = At{x — To). For simplicity, we will consider a spherical Gaussian profile 
given by 

t{x) = exp(—a;^/2i?^), x = |T|, (7) 

where R is the Gaussian width (in the case of point sources convolved with a Gaus¬ 
sian beam, R is the beam width). We could easily consider other functional profiles^ 
without any loss of generality. The expected number density of maxima per inter¬ 
vals (T, X + dx), {v, V -I- dv), ( k , k + dn) and (e, e -I- de), given a source of amplitude 
A in such spatial interval, is 


n{u, K,e\vs) 


8v^ 


rib 




1 /.. .. \2 a ,2 (.s-2res 


( 8 ) 


where i^s = A/ao is the normalized amplitude of the source, Ks = —AR'Ioi is the 
normalized curvature of the source and t" is the second derivative of the source 
profile T with respect to x at the position of the source. Note that in equation 
(8) we are taking into account that the shear of the source is zero since we are 
considering a spherical profile. It is useful to dehne a quantity ps that is related to 
the curvature of the source: 


Vs = 


P 


Ks = 


VsVs 


(9) 


^For example, more complicated beams or sources that are not pointlike but have some resolved 
structure. 


7 



3 THE DETECTION PROBLEM 


Equations (4) and (8) can be used to decide whether a source is present or not in a 
data set. The tool that allows us to decide whether a point source is present or not 
in the data is called a detector. In this section we will describe the Neyman Pearson 
detector (NPD). We will study its performance in terms of two quantities: the 
number of true detections and the number of false (spurious) detections that emerge 
from the detection process. Our approach fixes the number density of spurious 
detections and determines the number density of true detections in each case. 

3.1 The region of acceptance 

We consider a peak in the 2D data set characterized by the normalized amplitude, 
curvature and shear (i^, k, e). The number density of background maxima nh(^, k, e) 
represents the null hypothesis Hq that the peak is due to the background in the 
absence of source. Conversely, the local number density of maxima n(i^, k, e) rep¬ 
resents the alternative hypothesis, that the peak is due to the source added to the 
background. The local number density of maxima n{iy, k, e) can be calculated as 

COO 

n{v,K,e)= / dvsp(ya)n(y,n,e\vs)- ( 10 ) 

Jo 

In the last equation we have used the a priori probability p(vs), that gives the 
amplitude distribution of the sources. 

We can associate to any region TZ^:{v,K,,e) in the {v,K,e) parameter space two 
number densities nl and n* 


nl = 

dn dn de ni,(v, k, e). 

(11) 

Jn. 


n*= 1 

r 

dv dn den(v, k, e). 

(12) 


Jn, 

where nl is the expected number density of spurious sources, i.e. due to the back¬ 
ground, in the region 7?.* {n, k, e), whereas n* is the number density of maxima 
expected in the same region of the {n, k, e) space in the presence of a local source. 
The region 7?.* will be called the region of acceptance. 

In order to define the region of acceptance 7^* that gives the highest number den¬ 
sity of detections n* for a given number density of spurious detections nl , we assume 
a Neyman-Pearson Detector (NPD) using number densities instead of probabilities 

Tt ^ _ n{n,K,e) 

L{v, K, e) = — - - > L*, (13) 

nb(n, K, e) 

where L* is a constant. The proof follows the same approach as for the standard 
Neyman-Pearson detector. If L > L* we decide that the signal is present, whereas 



if L < L* we decide that the signal is absent. From this ratio L > L:^, we derive the 
region of acceptance, that is given by the sufficient linear detector ip (see Appendix) 

7^* : <p{v,k) > (/?*, (14) 

where is a constant and Lp{v, k) is given by 

(p{v,K) = av A-bn, a=- - 5 -, o=--(15) 

1 — 1 — 

We remark that the detector is independent of the shear e. This is due to the fact 
that we are considering a source with a spherical profile with shear Cg = 0. If the 
profile is not spherical, the detector may depend on the shear. 


3.2 Spurious sources and real detections 

Given a region of acceptance 7^*, we can calculate the number density of spurious 
sources and the number density of detections as given by equations ( 11 ) and ( 12 ) 


73. 


Ub 




dn (k^ — 1 + 


e )e 2 erfc (M). 
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(16) 
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7^ 


pOO pOO 

/ dvsp(vs) / (iK(7 — 1 + e“” )e“7'«-'2sy3) (Q) 

Jo Jo 

PVs - 1 


Q = M + Vg 


72(1-P")' 


(17) 

(18) 


Our approach is to fix the number density of spurious detections and then to deter¬ 
mine the region of acceptance that gives the maximum number of true detections. 
This can be done by inverting the equation (16) to obtain = ipt,{nl/nb; p,yg). 
Once ip* is known, we can calculate the number density of detections using equation 
(17). 


4 THE FILTERS 

Detection can, in principle, be performed on the raw data, but in most cases it is 
convenient to transform first the data in order to enhance the contrast between the 
distributions nb(v., k, e) and n{h', k, e). Hopefully, such an enhancement will help the 
detector to give better results (namely, a higher number of true detections). In this 
paper we will focus in the use of linear filters as a means to transform the data in 
such a way. Filters are suitable for this task because background fluctuations that 
have variation scales different from the source scale can be easily Hltered out while 
preserving the sources. Different Hlters will improve detection in different ways: this 
paper compares the performance of several filters. The filter that gives the highest 
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number density of detections, for a fixed number density of spurious sources, will 
be the preferred filter among the considered filters. 

Let us consider a filter 'i>{x;R,b), where R and b define a scaling and a translation 
respectively. Since the sources we are considering are spherically symmetric and we 
assume that the background is statistically homogeneous and isotropic, we will 
consider spherically symmetric filters, 


^'(f; R, b) = • (19) 

If we filter our background with 'I'(ai; R, 6), the filtered field is 

w{R,b) = j dx ^{x)'^{x] R,b). (20) 

The filter is normalized such that the amplitude of the source is the same after 
filtering: 

/dfr(f)«(f;R.0) = l. (21) 

For the filtered field equation (3) becomes 

pOO 

al = 27T dq q^+^^P{q)^\q). ( 22 ) 

Jo 

The values of p, 9m, 9c and all the derived quantities change accordingly. The 
curvature of the filtered source Kg can be obtained through equation (9), taking 
into account that for the filtered source 

pOO 

= dq q^T{q)il}{q). (23) 

JO 

Note that the function ■i/'(<z) will depend as well on the scaling R. As an application 
of the previous ideas, we study the detection of point sources characterised by a 
Gaussian profile t{x) = exp(—x^/2i?^), x = |x|, and Fourier transform T(q) = 
R? exp(—((7i?)^/2). This is the case we find in CMB experiments, where the profile 
of the point source is given by the instrumental beam, that can be approximated 
by a Gaussian. 

This profile introduces in a natural way the scale of the source R, the scale at 
which we filter. However, previous works in ID using the MHW, MF, SAF and the 
BSAF have shown that the use of a modified scale aR can significantly improve the 
number of detections (Gayon et al. 2000, Vielva et al. 2001a,b, Lopez-Ganiego et 
al. 2004a,b). Therefore, we generalise the functional form of these filters to 2D and 
allow for this additional degree of freedom a. 
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4.1 The matched filter (MF) 

Let us introduce a circularly-symmetric filter 5). The filtered field is given 

by equation (20). Now, we express the conditions to obtain a filter for the detection 
of the source s{x) = At{x) at the origin taking into account the fact that the source 
is characterized by a single scale Ro- We assume the following conditions: 

1. {w{Ro, 0)) = s(0) = A, i.e. w{Ro, 0) is an unbiased estimator of the amplitude 
of the source. 

2. The variance of w{R, b) has a minimum at the scale Ro, i.e. it is an efficient 
estimator. 

Then, the 2D filter satisfying these conditions is the so-called matched filter. As 
mentioned before, we will allow the hlter scale to be modified by a factor a. If 
a = 1 we have the well-known standard matched filter use in the literature. For a 
source with a Gaussian profile, a scale-free power spectrum P{q) oc q~'^ and allowing 
the filter scale to vary through the a. parameter, the modified matched filter is 

i’Mpiq) = N{a)z~^e~"^^ , z = qaR, (24) 


where 


m = 


2 - 1-7 

2 


N{a) 


1 1 ^ 2a^ 

A™7rr(TO)’ (1-I-a^) ’ 


(25) 


and r is the standard Gamma function. The parameters of the filtered background 
and source are 


p{a) = p = . , effia) = aR\ , ys{a) = pA (26) 

\ I + m \ I + m 

The corresponding threshold as compared to the standard matched filter (a = 1) 
is 

(27) 

^MF(a=l) 

where 

(28) 

We remark that for the standard matched filter the curvature does not affect 
the region of acceptance and the linear detector k) is reduced to (p = p. 


4.2 The scale-adaptive filter (SAF) 

The scale-adaptive filter (or optimal pseudo-filter) has been proposed by Sanz et al. 
(2001). The filter is obtained by imposing an additional condition to the conditions 
that define the MF: 

3. w{R, 0) has a maximum at {Ro, 0). 
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Assuming a scale-free power spectrum, P{q) (x q a modified scale aR and a 
Gaussian profile for the source, the functional form of the filter in 2D is 


i’SApiq) = N{a)z~*e 


2 t 2 

7 -I-z 

m 


z = qaR, 


(29) 


where 


N{a) 


I 1 
A™ 7rr(m) 7 -|- 


(30) 


and where m and A are defined as in equation (25), t is defined as in equation (28). 
The parameters of the filtered background and source are 


p{a) = p = 


m 


H, 


1 TTl -sj H2 

, ~m 

Vs{a) = 


, 0m{a) = aR 


2 Hi 
l + m\ H 3 ’ 


^2^7 + c(l -f m)A 


1 -I- TO V 7 -I- ctoA ’ 

where c = 2 t/m and 

Hi = 7^-I-270(1-I-to)-I- c^(l-I- to)(2-I-to), 

H2 = 7^-I-27CTO-I-c^to( 1-f to), 

H 3 = + 2 jc {2 + m) + c ^{2 + m){3 + m). 


(31) 

(32) 


(33) 


The corresponding threshold as compared to the standard matched filter (a = 1) 


IS 


^MF{a=l) 


a* 2 A™ ( 7 -I-CTO A) 

7^ ■ 


(34) 


4.3 The Mexican Hat wavelet (MH) 

The MH is defined to be proportional to the Laplacian of the Gaussian function in 
2D real space 

i’MHix) o: {1 — , a: = |a;|. (35) 

Thus, in Fourier space we get the modified Mexican Hat wavelet introducing the a 
parameter as follows 


-ipMHiq) = N{a)z^e z = qaR, 


(36) 


N{a) 


1 

TT 



Thus, the filtered background and source parameters are 


(37) 


p{a) = p = 



0m (a) = aR 



ys{a) 


2 

\/(2 + 0(3 + t ) 


A, 


(38) 
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where m and A are defined as in equation (25) and t is defined as in equation (28). 
The corresponding threshold is 


v{a) a* 

VMF{a=i) \/r(m)r(2 + 1) 


(39) 


4.4 The biparametric scale adaptive filter (BSAF) 

Lopez-Caniego et al. (2004b) have shown that removing condition 3 defining the 
SAF and introducing instead the condition 

3. w{Ro, b) has a maximum at {Ro, 0), 

leads to the new filter 

[1 + ciqR)"^] , (40) 

where c is an arbitrary constant related to the curvature of the maximum. For the 
case of a scale-free spectrum and allowing for a modified scale aR, the filter is given 
by the parameterized equation 


ipBSAF{q) = N{a)z^e (l-h cz^) , z = qaR, (41) 


N{a) 


I 1 1 

A™ TT F(m) 1 -I- ctoA ’ 


(42) 


where m and A are defined as in equation (25). We remark that c = 0 leads to the 
MF and if c = 2t/m'-f, with t defined as in equation (28), the BSAF defaults to the 
SAF. The parameters of the filtered background and source are 


p{a) = p = 


m 


Di 


1 m ^ D 2 


; ^m(o^) — CxR 


2 /i>r 
1 -I- m V D3’ 


ys{a) = 


1 -I- m V D 3 


ID 2 ^ 1 -I- c(l -I- to) A 


1 -I- ctoA 


where 

Di = 1-I-2c(l-I- to)-I- c^(l-I- to)(2-I-to), 

D 2 = 1 + 2cm + c^m{l + m), 

D3 = 1 + 2c{2 + m) + c^{2 + m){3 + m). 

The equivalent threshold is given by 

F{a) a*“^A’"(7 -I-ctoA) 

^MF{a=l) 


(43) 

(44) 


(45) 


(46) 
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5 ANALYTICAL RESULTS 


In this section we will compare the performance of the different filters previously 
introduced. We use as example the interesting case of white noise as background. 
This is a fair approximation to the case presented in Figure 1, where the sources 
are embedded in a background that is a combination of instrumental noise (approx¬ 
imately Gaussian) and a small contribution of residual foregrounds that have not 
been perfectly separated. For this example, we will consider sources with intensities 
distributed uniformly between zero and a upper cut-off. 

The comparison of the filters is performed as follows: we fix the number density 
of spurious detections, the same for all the filters. Then, for any given filter we 
calculate the quantities (t„, p and j/g. Using equation (16) it is possible to calculate 
the value of lyj* that defines the region of acceptance. Then we calculate the number 
density of real detections using equation (17). The filter that leads to the highest 
number density of detections will be the preferred one. We do this for different values 
of a in order to test how the variation of the filtering scale affects the number of 
detections. 

5.1 A priori probability distribution 

As mentioned before, we will test a pdf of source intensities that is uniform in the 
interval 0 < A < Ac- In terms of normalized intensities, we have the pdf 

P{’^s) = —, (47) 

We will consider a cut-off in the amplitude of the sources such that t'c = 2 after 
filtering with he standard MF, that is, we will focus on the case of faint sources 
that would be very difficult to detect if no filtering was applied. Note that while the 
value Vc is different for each filter (because each filter leads to a different dispersion 
(To of the filtered field), the distribution in source intensities A is the same for all 
the cases. 

5.2 Results for white noise 

We want to find the optimal filter in the sense of maximum number of detections. 
For the sources, we use a uniform distribution with amplitudes in the interval A G 
[0, 2]cto, where ctq is the dispersion of the linearly-filtered map with the standard 
MF. We focus on the interesting case of white noise (7 = 0) and explore different 
values of nj and R. 

We study the performance of the different filters as a function of a. This allows 
us to test how the variation of the natural scale of the filters helps the detection. In 
the case of the BSAF, which has an additional free parameter, c in equation (41), 
for each value of a we determine numerically the value of c that gives the highest 


14 




Figure 2: The expected number density of detections n* as a function of a for 7 = 0 for 
the BSAF (c has been obtained by maximising the number of detections for each value 
of a), MF and SAF filters. We consider the case R = 1.5, = 0.01 






Figure 3: The expected number density of detections n* as a function of a for 7 = 0 for 
the BSAF (c has been obtained by maximising the number of detections for each value 
of a), MF and SAF filters. We consider the case R = 2, nl = 0.01 





R 


a 

c 

''^BSAF 

n*MF 

RD[%] 

1.5 

0.005 

0.5 

-0.44 

0.0507 

0.0484 

4.7 


0.01 

0.5 

-0.46 

0.0709 

0.0620 

14.3 

2 

0.005 

0.4 

-0.54 

0.0396 

0.0335 

18.2 


0.01 

0.4 

-0.54 

0.0567 

0.0406 

39.6 

2.5 

0.005 

0.3 

-0.64 

0.0320 

0.0245 

33.3 


Table 1: Number density of detections n* for the BSAF and the standard MF(q; = 1) 
with optimal values of c and a for different values of and R. RD means relative 
difference in number densities in percentage; RD = 100(—1 + n*^sAFl'^*MF)- 


number of detections. Then the BSAF with such c parameter (that is a function of 
a, nl and R) is compared with the other filters. 

In Figure 2, we plot the expected number density of detections n* for different 
values of a, R = 1.5 pixels and nl = 0.01. Note that for the 2D case the MHW 
and SAF are the same filter for 7 = 0 and we have only included the latter in our 
figures. In this case, the curve for the BSAF always goes above the other filters. 
The maximum number of detections is found for small values of a. In this region, 
the improvement of the BSAF with respect to the standard matched filter is of the 
order ~ 15%. 

In Figure 3, we show the results for i? = 2. We have increased the beam width 
as compared with the previous example and leave unchanged the number density of 
false detections. The BSAF outperforms all the other filters, and for small values 
of a the improvement is of the order ~ 40%. Note that in this figure the MF takes 
values a e [0,1]. For greater values of a, with R = 2 and nl = 0.01, we can not 
solve for in the implicit equation (16) and can not calculate n*. 

We remark that filtering at scales much smaller than the scale of the pixel does 
not make sense. This is due to the fact that we are not including the effect of the 
pixel in our theoretical calculations and, thus, the results would not exactly follow 
what would be found in a real image. Therefore, we present the results only for 
those values of a such that aR is at least ~ 1. 


6 CONCLUSIONS 

Several techniques have been introduced in the literature to detect point sources in 
two-dimensional images. Examples of point sources in astronomy are far galaxies as 
detected by CMB experiments. An approach that has been thoroughly used in the 
literature for this case consists in linear filtering the data and applying detectors 
based on thresholding. Such approach uses only information on the amplitude of the 
sources: the potentially useful information contained in the local spatial structure 
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of the peaks is not used at all. In our work, we design a detector based on peak 
statistics that uses the information contained in the amplitude, curvature and shear 
of the maxima. These quantities describe the local properties of the maxima and 
are used to distinguish statistically between peaks due to background fluctuations 
and peaks due to the presence of a source. 

We derive a Neyman-Pearson detector (NPD) that considers number densities 
of peaks which leads to a sufficient detector that, in the case of the spherically 
symmetric sources that we consider, is linear in the amplitude and curvature of the 
sources. For this particular case, then, the information of the shear of the peaks is 
not relevant. In other cases, however, it could be useful. 

It is a common practice in astronomy to linear filter the images in order to 
enhance very faint point sources and help the detection. The best filter would be 
the one that makes easier to distinguish between peaks coming from the background 
alone and those due to the presence of a source, according to the information used by 
the detector. In the case of simple thresholding, that considers only the amplitude 
of the peaks, the answer to the question of which is the best filter (in the previous 
sense) is well known: the standard matched filter. But in the case of the Neyman- 
Pearson detector, that considers other things apart from mere amplitudes, this is 
not longer true. 

We have compared three commonly used filters in the literature in order to assess 
which one of them performs better when detecting sources with our scheme. In 
addition, we have designed a filter such that optimizes the number of true detections 
for a fixed number of spurious sources. The optimization of the number of true 
detections is performed by using the a priori pdf of the amplitudes of the sources. 
This filter depends on two free parameters and it is therefore called biparametric 
scale adaptive filter (BSAF). By construction, the functional form of the BSAF 
includes the standard MF as a particular case and its performance in terms of 
number of true detections for a fixed number of spurious detections must be at least 
as good as the standard MF’s one. 

Following the work done in the ID case, we generalize the functional form of 
the filters to 2D and introduce an extra degree of freedom a that will allow us to 
filter at different scales ai?, where R is the scale of the source. This significantly 
improves the results. 

We have considered an interesting case, a uniform distribution of weak sources 
with amplitudes A G [0, 2]cro, where uo is the dispersion of the field filtered with the 
standard matched filter, embedded in white noise (7 = 0 ). We have tested different 
values of the source size R and of the number density of spurious detections that 
we fix. We find that the BSAF improves the number density of detections up to 
~ 40% with respect to the standard MF (a = 1) for certain cases. Note that since 
the Neyman-Pearson detector for the standard MF (a = 1) defaults to the classic 
thresholding detector that is commonly used in astronomy, the results of this work 
imply that it is possible, under certain circumstances, to detect more point sources 
than in the classical approach. 

The generalization of these ideas to other source profiles and non-Gaussian back- 
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grounds is relevant and will be discussed in a future work. 
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A Appendix 

We will show in this Appendix that k) > given in equation (14) is a sufficient 
linear detector, i.e., the detector is linear in the threshold v and the curvature k and 
the data it uses is a sufficient statistic to decide if a peak is a source (independent 
of the a priori probability P(i^s)). The ratio L{iy,K,e\vs) = n{v,K,e\vs)/'nb{v,K,,e) 
can be explicitly written as 

The criterion for detection can be written as 

poo 

L{iy,K)= / (49) 

^0 

where L* is a constant. L is a function of (/?, 

/- \ _ I J, 1 ~ PVs Vs ~ P /'cn^ 

ip[v, k) = av + OK, a = -- b = --(50) 

I — I — 

By differentiating L with respect to tp we find that 

> 0, (51) 

and therefore setting a threshold in L is equivalent to setting a threshold in tp: 

L{iy, k)> ‘p{v, k) > ipt, (52) 

where k) is given by equation (50) and (p* is a constant. 


dp Jo 
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